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Abstract 

The Sweet-Parker layer in a system that exceeds a critical value of the Lundquist number (5) 
is unstable to the plasmoid instability. In this paper, a numerical scaling study has been done 
with an island coalescing system driven by a low level of random noise. In the early stage, a 
primary Sweet-Parker layer forms between the two coalescing islands. The primary Sweet-Parker 
layer breaks into multiple plasmoids and even thinner current sheets through multiple levels of 
cascading if the Lundquist number is greater than a critical value S'c ~ 4 x 10^. As a result 
of the plasmoid instability, the system realizes a fast nonlinear reconnection rate that is nearly 
independent of S, and is only weakly dependent on the level of noise. The number of plasmoids in 
the linear regime is found to scales as S^^^ , as predicted by an earlier asymptotic analysis (Loureiro 
et al, Phys. Plasmas 14, 100703 (2007)). In the nonlinear regime, the number of plasmoids follows 
a steeper scaling, and is proportional to S. The thickness and length of current sheets are found 
to scale as S~^, and the local current densities of current sheets scale as Heuristic arguments 
are given in support of theses scaling relations. 
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I. INTRODUCTION 



Recent studies of nonlinear reconnection in large high-Lundquist-number (S) plasmas, 
based on resistive magnetohydrodynamics (MHD) [l| as well as fully kinetic simulations 
that include a collision operator have produced a surprise. It is seen in these studies that 
above a critical value of the Lundquist number, the system deviates qualitatively from the 
predictions of Sweet-Parker theory |3|, IJ] which has been the standard model for reconnection 
in the high-S* regime. In the Sweet-Parker model, the reconnection layer has the structure 
of Y-points, with a length of the order of the system size, and a width given hj 6sp = 
L/S^^'^, where S = LVa/t] is the Lundquist number based on the system size L, the Alfven 
speed Va, and the magnetic diffusivity 77. The Sweet-Parker model is usually considered 
to be a model of "slow" reconnection as it predicts the reconnection rate to scale as 5*"^^^. 
(The Petschek model jsl predicts a much weaker dependence on S, with the maximum 
reconnection rate ~ 1/ log S. However, it has become clear over the years that Petschek 
model is realizable only when the resistivity is locally enhanced around the reconnection 
site. gI) For weakly collisional systems such as the solar corona, the Lundquist number is 
typically very large (~ 10^^ — 10^^) and the Sweet-Parker reconnection time scale is of the 
order of years, much too slow to account for fast events such as solar flares. The Sweet- 
Parker model is based on the assumption of the existence of a long thin current layer. 
Although it has been known for some time that such a thin current layer may be unstable 
to a secondary tearing instability (referred to hereafter as the plasmoid instability) which 
generates plasmoids, 7, y, o], 10, [ll|, 1^ it has been realized only fairly recently that the 
Sweet-Parker layer actually becomes more unstable as the Lundquist number increases, with 
a linear growth rate 7 ~ 5*^/^ and the number of plasmoids ~ S^^^. l|, In a recent paper 
(Ref. hereafter referred to as Paper I), Bhattacharjee et al. have presented numerical 
results that suggest strongly that as a consequence of the plasmoid instability, the system 
evolves into a nonlinear regime in which the reconnection rate becomes weakly dependent 
on S. 

The primary goal of this paper is to strengthen the results obtained in Paper I in two 
significant ways: first, to present new simulation results with a modified initial condition 
that enables us to obtain stronger scaling results on the nonlinear reconnection rate and 
the number of plasmoids generated in the nonlinear regime, and second, a simple heuristic 
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model that is consistent with the results of the simulation and fortifies the claim in Paper 
I that the reconnection rate in the nonlinear regime of the plasmoid instability is fast and 
independent of S. 



II. NUMERICAL MODEL 

The initial condition in Paper I does not have a thin current sheet to begin with. It 
has four magnetic islands and is unstable to an ideal coalescence instability. After the 
onset of the coalescence instability, a Sweet-Parker current sheet is created when two islands 



are attracted toward each other. IJ] In this case, there is a relatively long initial transient 
period before the reconnection process starts. Furthermore, the dynamics of the system are 
complicated by the sloshing of coalescing islands that causes the primary Sweet-Parker layer 
to lengthen or shorten from time to time. This makes it difficult to verify the predictions 
of linear theory in this particular system. The present study seeks remedies to these two 
drawbacks. We still consider the merging of two islands, but now put them in close contact 
initially. There is an initial current layer between the fiux tubes, which quickly (typically 
within less than one Alfven time) adjusts its width depending on the Lundquist number 
to form the primary Sweet-Parker layer, which may subsequently become unstable to the 
plasmoid instability. In this new system, the transient period is significantly shortened and 
the sloshing between islands is largely eliminated. It is still not easy to measure the linear 
growth rate in this new model, but we can at least verify the scaling of the number of 
plasmoids in the linear regime. 

The initial condition is similar to the model of Uzdensky and Kulsrud. [l^ To start with, 
let tpQ = cos (ttx) sin (27r |2;|) /27r and Bq = V^/'o x y in the domain {x,z) G [—1/2,1/2] x 
[—1/2,1/2]. The ipQ so defined satisfies V^ipo = — Stt^-j/'o- If the pressure is set to po = 
C + 57r^^/;Q/2, with C an arbitrary constant, the system is in force balance. However, the 
magnetic field defined by tpo has a tangential discontinuity at 2; = 0, which causes numerical 
difficulties. We smooth it out as ip = tanh (a^) cos (ttx) sin (27r2;) /27r, where a is a large 
number. This smoothed function no longer satisfies V'^tp = f{ip) and the magnetic force 
cannot be balanced entirely by pressure. However, the magnetic force can be canceled to 
a large extent if the pressure is set to p = po + {B^q — -B^)/2, where B^o = —dzipQ and 
B^ = —dzipi respectively. We assume the isothermal equation of state, p = 2pT, in our 
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simulation. We choose T = 3, C = 2T, and a = 100. For these parameters, the initial 
plasma density varies from 0.96 to 1.1 and the plasma beta {/3 = p/B'^) obeys the inequality 
/3 > 6. The system is, therefore, nearly incompressible. Figure [T] shows the initial current 
density and magnetic field lines. 

We find in the present study that the plasmoid instability depends on the noise level of the 
system, at least when the Lundquist number is not far above the critical value. Due to the 
outflow in the primary Sweet-Parker layer, if the noise level is low, the plasmoid instability 
may not grow to visible size before being convected out. When we seed the system initially 
with random noise, at relatively low values of the Lundquist number (~ 10^) we obtain a 
short burst of plasmoids, following which the current layer becomes stable again after all 
the plasmoids are convected out. Because of this, we have included a random forcing in the 
system, which enables us also to study the effect of noise level on the reconnection rate. 
(Notice that we did not apply any random forcing or noise in Paper I. The sloshing between 
coalescing islands is a natural source of noise, but uncontrolled.) The governing equations 
for the time evolution of the system are: 

9iP + V-(pv) = 0, (1) 

5t(pv) + V ■ (pvv) = -Vp - VijV^tp + ef (x, t), (2) 

atV' + V ■ V^/' = r/V^V^, (3) 

where a random forcing term ef(x, t) is added to the right hand side of the momentum 
equation The forcing function is white noise in both space and time with (f) = 0, 
and e is a small parameter which controls the noise level. By white noise we mean that 
(/j(x, t)/,(x', t')) ~ 5ij5(x — x.')S(t — t'), where i and j can be x or 2 and ( ) is the ensemble 
average. Care has to be taken to ensure that the discrete representation is independent of 



the time step. Our implementation is similar to that in Ref. IGj. It is convenient to set 
f = pa, where a is a random acceleration. In a single time step, the momentum density 
evolves from pv to pv + epaAt (neglecting other forces), and the kinetic energy density 
evolves from pv'^/2 to pf ^/2 + epa ■ vAt + e^pa^At^/2. Therefore, the average power density 
from the random force is 

(.pa.vAt + eVW/2)^l , 
At 2 ^ \ 
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where (a ■ v) = is used, and = a ■ a. At each grid point, we set a^, az to random 
numbers between —1 and 1 with a uniform probabihty distribution, divided by v^At. That 
is, 

rand(— 1, 1) 
a,; = - 



'At 

Then At {a?) = 2/3 and is independent of At. The average power density is and the 

total power e^M/3, with M the total mass (M ^ 1 in our simulation). We use e = 10~^— 10~^ 
in our simulations and the corresponding power density ranges from 3 x 10"^^ to 3 x 10~^. 
Since our simulations typically last only a few Alfven times, the energy input from random 
forcing is negligible compared to the total magnetic energy (~ 0(1)) in the system. This 
ensures that the random forcing only provides noise for the instability to grow but does not 
otherwise alter the system in a significant way (see more discussion in Sec. Illll) . 

Our numerical algorithm [l^] uses finite differences with a five-point stencil in each di- 
rection, and a second-order accurate trapezoidal leapfrog method for time stepping. We use 
a uniform grid along the x direction, and a nonuniform grid in the z direction that packs 
high resolution around z = in order to resolve the sharp spatial gradients in the recon- 
nection layer. Perfectly conducting [dtil) = 0), impenetrable (v • n = 0), and free slipping 
(n ■ V(n X v) = 0) boundary conditions are assumed (n is the unit normal vector to the 
boundary). We further assume reflection symmetry along the x axis and only the region 
2 > is simulated. The highest resolution is 12800 in x and 1600 in z, with the smallest 
grid size Az = 5.4 x 10~^. 



III. NUMERICAL RESULTS 



One of the key objectives of this study is to determine the scaling of reconnection rate 
in the plasmoid- unstable regime. To quantify the speed of reconnection, we measure the 
time it takes to reconnect a significant portion of the magnetic flux within the two merging 
islands. The amount of magnetic flux in an island is ipmax —"^s, where ipmax is the maximum 
of ip in the island and ips is the value of at the separatrix separating the two merging 
islands. Initially ipmax — 0.16 and it remains approximately unchanged since the resistivity 
is low; therefore it suffices to just measure ips- We denote the time it takes to reconnect from 
ips = 0.01 to tps = 0.05 as tree- The starting point tps = 0.01 is chosen to allow the plasmoid 
instability to build up, and the end point ips = 0.05 is chosen such that the reconnection 
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layer does not shorten too much compared with that in the initial condition. The range 
corresponds to reconnecting 25% of the initial flux. 

Figure [2] shows the reconnection time t^ec for various S and e. For lower S, the reconnec- 



tion time tree ~ 5*^/^, as expected from the Sweet-Parker theory. ISj The critical Lundquist 
number for plasmoid instability is about S'c — 4 x 10^. Above Se, the reconnection time 
tree deviates from the Sweet-Parker scaling and becomes shorter. In the plasmoid unstable 
regime, the reconnection time is nearly independent of S. However, the reconnection time 
has a weak dependence on the noise level throughout the S range we have tested. The 
plateaued values of tree in the high-S' regimes for e = 10~^, 10~^, and 10~^ are 5.30 ± 0.27, 
6.10 ± 0.41, and 7.05 ± 0.16, respectively. Here we take the average values over the range 
5 = 5 X 10^ to 3 X 10^, and the errors represent the standard deviation. We have tested 
the convergence of our numerical results by varying the resolution, the time step, and the 
random seed for selected runs. These are represented by multiple data points with the same 
parameters in Figure O The results are fairly consistent, with fluctuations no more than a 
few percent. The dependence of tree on e may be tentatively fit with a power law, which gives 
tree ~ £-0 062 Howcvcr, glvcu the limited range of the parameter space we have explored, 
this scaling should not be considered as conclusive. 

The global characteristic values for Va and B are about 1, which yield the normalized 
average reconnection rate as 

1 /#A 0.04 



BVa \ dt I tr 

In the high Lundquist number regime, tree — 5 to 7 from our simulations and the normalized 
reconnection rate is in the range 0.006 to 0.008. The normalized reconnection rate obtained 
here is similar to the result of Paper I. 

Figure [3] shows a time sequence of the current density, overlaid with magnetic field lines, 
within a small area (1/1000) of the whole domain for a case with e = 10~^ and S* = 3 x 10^. 
The initial current layer (panel (a)) quickly thins down to form the primary Sweet-Parker 
layer (panel (b)), which becomes unstable to the plasmoid instability (panel (c)). As the 
instability proceeds, the plasmoids grow in size and the current sheets between plasmoids 
are again Sweet-Parker like (panel (d)). These secondary Sweet-Parker current sheets are 
thinner than the primary one and are again unstable to the tertiary plasmoid instability 



(panel {e_ 



in Ref. 



I. This process of multiple stages of cascading resembles the scenario envisaged 



ll| . The plasmoids can merge to form larger ones and new plasmoids are constantly 



generated (panel (f)). Figure H] shows the global configuration at a later time t = 3.9, as 
well as close-ups of the reconnection layer. The figure shows that on a large, coarse-grained 
scale, the configuration looks Sweet-Parker like, except for the important difference that the 
reconnection layer is no longer a single extended current sheet, but is made up of a sequence 
of copious plasmoids and current sheets. 

Linear theory predicts that the number of plasmoids, n^, scales as S"^/^.!!!, [l^ We verify 
this by counting the maximum number of plasmoids within the central region, —0.25 < x < 
0.25, before the plasmoid instability becomes highly nonlinear (roughly corresponds to panel 
(c) in Figure [3]). Figure O shows the number of plasmoids versus S, for e = 10~^, in both 
linear and nonlinear (see the discussion later) regimes. The result in the linear regime is 
in good agreement with the S^^^ scaling predicted by asymptotic analysis. This scaling has 
been verified by Samtaney et al. with local simulations up to 5 = 10*. j 



In the fully nonlinear regime, the plasmoid dynamics are very complicated and constantly 
evolving. Plasmoids may grow in size, coalesce with each other to form larger plasmoids, 
and finally get ejected into the downstream region. Meanwhile, new plasmoids are con- 
stantly generated in the reconnection layer. We may regard the reconnection layer with 
multiple plasmoids as a statistical steady state. As a simple, first approximation, we expect 
the cascading to stop when the current sheet segments between plasmoids become stable. 
We may imagine the reconnection layer as a chain of plasmoids connected by marginally 
stable Sweet-Parker current sheets. For given rj and Va, the critical length of a marginally 
stable current layer is = ScTj/Va- Therefore we expect the number of plasmoids in the 
nonlinear regime, Up^, to scale like L/Lc ~ S/Sc- Furthermore, the thickness of each 



Sweet-Parker sheet is 6c ~ Lc/ ~ vVS^/^a ~ ^spy Sc/ S , and the current density 
J B / 5c ^ BVA/Tj^fSl ~ BS/L\fS'c ■ If we identify the reconnection rate with 77 J, then 
the reconnection rate ~ rjj ~ BVa/ \fS'ci which is independent of S. This is consistent 
with our finding that the reconnection rate is nearly independent of S in the high-S" regime. 
Clearly, the assumption that all current sheets are marginally stable, and therefore all iden- 
tical, is oversimplified. If we look at the individual current sheets, there are a whole variety 
of them, each with a different length, width, and current density. Therefore, the system is 
better described with a statistical approach. If we neglect complications such as asymmetry 
or background shear flow, and consider the simple Sweet-Parker picture for each current 
sheet, then the local Lundquist number Siocai = VaI/t] is the only dimensionless parameter 
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associated with it. Here / denotes the length of the current sheet. The current sheet thick- 



ness will be 5 ~ 1/ VSiocah as predicted by the Sweet-Parker theory. The local Lundquist 
number being greater or smaller than Sc determines whether a current sheet may or may 
not further break into plasmoids and even smaller current sheets. Because it is the local 
Lundquist number that determines the cascading of a local current sheet to even smaller 
scales, we hypothesize that the probability distribution of Siocai is independent of the global 
Lundquist number. The underlying assumption is that, if we consider the ensemble of local 
current sheets and characterize each current sheet by a dimensionless parameter Siocau there 
is a similarity across systems of different global Lundquist numbers. If we further assume 
that the local upstream Alfven speed is determined by global conditions, it follows that 
statistically the length / and the thickness 5 of a current sheet scale as rj , and the current 
density J scales as ri~^ . If we consider the simple picture that two neighboring current sheets 
are separated by a plasmoid, then / ~ rj implies the number of plasmoids in the nonlinear 
regime ri~^ . 

Now we proceed to examine whether the conclusions from the simple heuristic argument 
are consistent with our simulation data. Here we present the results from cases with e = 10^'^. 
Results from other values of e are similar. We count the number of plasmoids by first 
identifying X-points and 0-points along 2 = 0. There are two types of 0-points, the local 
minimum (type I) and the local maximum (type II) of ip. Likewise, there are two types of 
X-points, one with dlip < 0, dlip > (type I) and the other with d^^jj > 0, d'^tp < (type 
II). In the linear regime, only X-points and 0-points of type I are present. When plasmoids 
start to coalesce with each other, type II null points may be created. We count the number 
of type I 0-points within —0.25 < x < 0.25 as the number of plasmoids in the nonlinear 
regime. As shown in Figure the number of plasmoids in the nonlinear regime appears to 
agree with the ~ S scaling. Because the number of plasmoids fluctuates, the median value 
is used; the error bar indicates the range between the first and the third quartiles. Notice 
that although the number of plasmoids in the nonlinear regime follows a steeper scaling than 
that in the linear regime, it is not until about S = 2 x 10^ that the nonlinear scaling catches 
up with the linear counterpart. This is because at lower S, the coalescence and ejection of 
plasmoids exceeds the generation of new plasmoids. Equating the estimate, ^ S'^^^ /27r, 



from the linear theory l|, ll3| with the heuristic nonlinear estimate, ~ S/ Sc, and using 
Sc — 4: X 10^, we obtain ^ rip^ when S* ~ 1.2 x 10^, which is in approximate agreement 
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with the observed 5* = 2 x 10^. 

To examine the statistics of current sheets, we have to first set up a diagnostic for a 
current sheet, which is subject to a certain degree of arbitrariness. We search for local 
maxima of J greater than 10% of the global maximum Jmax within —0.25 < x < 0.25 as 
potential sites of current sheets. However, two neighboring maxima are regarded as separate 
current sheets only when the trough between them is lower than 25% of the greater of the 
two. The length I and the thickness 5 of a current sheet are measured by the locations where 
the current density drops to 25% of the local maximum J of the current sheet. 

Figures E] and [7] show scalings of the thickness 5, length / , and current density J with 
respect to the global Lundquist number S. The data are collected from time slices during the 
period to reconnect 25% of the flux in each case. Again the median values are used, and error 
bars indicate the range between the first and the third quartiles. Also shown for reference are 
the predictions from the heuristic argument based on marginally stable current sheets, i.e., 
/ ~ ~ LSc/S, 5 ~ 5sp{Sc/Sf'^ ~ LSl'^/S, and J ^ BS/LSl'^. It is evident that 

the characteristics of current sheets are distributed over a broad range, as indicated by the 
rather large error bars. Clearly, the observed quantities follow the expected scalings. Quite 
surprisingly, the predictions for 5 and J from the heuristic argument are in good agreement 
with the observed median values, even though the argument itself is rather crude. However, 
the prediction of I appears to be systematically an overestimate, and lies at the larger end 
of the numerically observed lengths. This is consistent with the fact that Lc = LSc/S is 
the critical length just above which the plasmoid instability is triggered. One may wonder 
how the prediction of 6 can be consistent with the observed values when the prediction 
of / is an overestimate. A possible explanation is that, the heuristic argument assumes 
Sweet-Parker-like local current sheets, but clearly not all current sheets in simulations are 
Sweet-Parker-like. The existence of non-Sweet-Parker-like current sheets is evident from the 
movie available online, which is for the case S' = 3xl0^,e = 10"^. 

Let us now take a more detailed look into the statistics of current sheets. Figure [H] shows 
the probability distribution functions (PDFs) of r]J for S = 10^, 2 x 10^, 3 x 10^, from 
cases with e = 10"'^. The case S = 10^ has been done with two runs. The PDFs of rjJ from 
different runs clearly show a degree of similarity, which lends some support to our hypothesis 
of similarity across systems of different global Lundquist number. However, we also notice 
some differences between the PDFs from different runs. Even the two runs with S = 10^ 
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show a significant variation in the PDFs. Therefore, more study is needed to further assess 
the vahdity of our hypothesis. Ideally the same global setting should be repeated many times 
with different random seeds for better statistics, but that is computationally too expensive 
to be done at the present time. 

Before we conclude this Section, we remark on a few subtle issues. In the heuristic 
argument given above, we have used the quantity rjj to estimate the reconnection rate. 
Strictly speaking, this is valid when the X-point and the stagnation point of the flow coincide, 
which is not necessarily the case when the reconnection layer is embedded with multiple X- 
points and plasmoids. Notwithstanding this caveat, we generally find that the peak value 
of 7] J is a reasonable measure of the reconnection rate. 

We also address the issue of whether random forcing, by itself, can significantly enhance 
the reconnection rate. An estimate of the effect of reconnection rate due to random fluc- 
tuations is the quantity |v x B| at the reconnection layer, where v is the random velocity 
fluctuation. Our estimates indicate that the contribution of random fluctuations is less than 
1% of the observed reconnection rate in the high-S* regime. This conclusion is reinforced 
by the fact that in the plasmoid stable regime, the variations in tree for different e are neg- 
ligible. This is qualitatively different from a recent turbulent magnetic reconnection study 
by Loureiro et al, 19| where the system is more strongly driven, and the reconnection rate 
shows a noticeable dependence on the magnitude of the forcing even for a Lundquist number 
as low as 10^. 



IV. SUMMARY AND CONCLUSION 

In summary, we have shown through a series of simulations that resistive MHD can 
achieve a fast reconnection rate in the high-Lundquist-number regime. Fast reconnection is 
facilitated by the plasmoid instability. The resultant reconnection rate is independent of 5* 
and is weakly dependent on the noise level. We have verified the 5''^/^ scaling of the number 
of plasmoids in the linear regime, as predicted in Refs. In the nonlinear regime, the 

number of plasmoids follows a steeper scaling and is proportional to 5*. We also have done 
statistical studies of the local current sheets, and found that the current sheet thickness and 
length both scale as S~^, while the current density scales as S. These findings are consistent 
with our heuristic argument and the claim that the reconnection rate is independent of S 
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in the high-S* regime. 

The fast reconnection rate we have obtained is approximately O.OIVa-B, which is similar 
to the values from other recent resistive MHD studies, [2^ ^lut is smaller than the typical 
reconnection rate from collisionless two-fluid or particle-in-cell simulations by an order of 
magnitude. Which rate will be realized depends on how collisional the system is. If the 
Sweet-Parker thickness Ssp is greater than the ion skin depth di (or the ion Larmor radius 
at the sound speed, ps, if there is a guide field) in a system, a Sweet-Parker layer is likely 
to form first. On the other hand, if 6 qp < dj (or Oo), the reconnection would likely proceed 



dominated by collisionless effects. 
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23|, 
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26| An interesting regime that has not 



drawn much attention before is when 5sp > di (or p^) but S > Sc- Then we expect the 
plasmoid instability to set in and the primary Sweet-Parker layer will break into segments. 
This brings the thickness further down to 5 ~ 6sp{Sc/SY^'^. li 5 > di (or ps) then the system 
is still dominated by collisional effects and we may end up getting a reconnection rate of 
O.OIVa-B. However, ii 6 < di (or ps) then the system will be in the collisionless regime. It is 
possible the reconnection rate will be further enhanced. Indeed, the recent particle-in-cell 

CI 

simulation by Daughton et al. suggests this possibility [2| but more needs to be done to 
determine if this is a general trend. 
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Figure 2: (Color online) The reconnection time tree for various S and e. The dashed line is the 
Sweet-Parker scaling. 
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Figure 3: (Color online) Time sequence of the current density, overlaid with magnetic field lines, 
for the case 5 = 3 x 10^, e = 10^'^. The dashed line indicates the separatrix. The dots and crosses 
indicate O points and X points, respectively. The blue and the red colors indicate the two types 
of X and O points (blue for type I and red for type II) . Notice that the vertical axis is stretched 
for better visualization. 
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t= 3.90e+00, J [-2.646+04. 1.40e+04] 
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Figure 4: (Color online) The global configuration at a later time, t = 3.9. Show on the right are 
close-ups of the reconnection layer. 



10' 



in 

(M 



(M 

o 

I 



10 

o 

E 
w 

Q. 

6 



^3/8 



O Linear 
• Nonlinear 



/ 

/ o 



rA- 



10" 



10 



10 

S 



10 



Figure 5: (Color online) Number of plasmoids in the linear and the nonlinear regime. 
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Figure 6: (Color online) Scalings of the current sheet thickness S and length I with respect to 
the global Lundquist number S. Also shown for reference are the predictions from the heuristic 
argument. 
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Figure 7: (Color online) The scaling of the local peak current density. Also shown for reference is 
the prediction from the heuristic argument. 
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Figure 8: (Color online) The probability distribution function of local rjj for different S", from cases 
with e = 10"^. 
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